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Abstract 

We propose a new approach for studying the notion of the instantaneous frequency of a signal. We 
build on ideas from the Synchrosqueezing theory of Daubechies, Lu and Wu and consider a variant of 
Synchrosqueezing, based on the short-time Fourier transform, to precisely define the instantaneous fre- 
quencies of a multi-component AM-FM signal. We describe an algorithm to recover these instantaneous 
frequencies from the uniform or nonuniform samples of the signal and show that our method is robust 
to noise. We also consider an alternative approach based on the conventional, Hilbert transform-based 
notion of instantaneous frequency to compare to our new method. We use these methods on several 
test cases and apply our results to a signal analysis problem in electrocardiography. 
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diography 
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1 Introduction 

In a recent paper [5], Daubechies, Lu and Wu proposed and analyzed an adaptive wavelet-based signal 
analysis method that they called "Synchrosqueezing." The authors considered continuous-domain signals 
that are a superposition of a finite number of approximately harmonic components, and they showed that 
Synchrosqueezing is able to decompose an arbitrary signal of this type. They showed that for such signals, 
Synchrosqueezing provides accurate instantaneous frequency information about their constituent compo- 
nents. 

The concept of instantaneous frequency is a natural extension of the usual Fourier frequency that describes 
how fast a signal oscillates locally at a given point in time, or more generally, the different rates of oscilla- 
tion at a given time. However, instantaneous frequency has so far remained a somewhat heuristic concept 
and has lacked a definition that is both mathematically rigorous and entirely satisfactory [TT] . The analysis 
of signals from samples spaced nonuniformly in time is also an important problem in several applications, 
arising in radar detection, audio processing, seismology and many other fields pi I12|. In this paper, we 
propose a new approach based on Synchrosqueezing to precisely define the instantaneous frequencies of a 
signal, and to recover these instantaneous frequencies from the uniform or nonuniform samples of the signal. 
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We consider a class of multi-component AM-FM signals, comparable to the type studied in [S], and define a 
variant of Synchrosqueezing based on ideas similar to those in [5] but using the short-time Fourier transform 
(STFT). We show that by applying this modified Synchrosqueezing transform to an impulse train weighted 
by the samples, we can determine the instantaneous frequencies of the signal with high accuracy. We show 
furthermore that this procedure is robust with respect to noise. For purposes of comparison, we also 
consider a parallel approach based on the conventional, Hilbert transform-based notion of instantaneous 
frequency, combined with a well known least-squares method for bandlimited signal reconstruction from 
nonuniform samples |10| . We apply both methods to several test cases and compare their performance. We 
also consider a problem concerning the extraction of respiration data from a single-lead electrocardiogram 
(ECG) signal, and show how these methods can be used to study it. 

This paper is organized as follows. In Section [21 we briefly review the conventional approach to instanta- 
neous frequency. In Section [3l we describe our Synchrosqueezing-based algorithm and present the main 
new concepts and results of the paper. Our adaptation of the least-squares method for bandlimited func- 
tions is described in Section 31 We then perform our numerical experiments in Section |S1 and discuss the 
application to ECG analysis. 



2 Background Material 

We first discuss the precise meaning of the instantaneous frequency (IF) of a signal and some of the obstacles 
encountered in computing it. We start by making a few definitions. Let / be a tempered distribution. We 
denote the forward and inverse Fourier transforms of / by / and /, using the normalization e^'^^^ — e~^^ . 
For a fixed window function g in the Schwartz class S, we denote the modified short-time Fourier transform 
of / by 

/OO 
/(x)g(x ~ t)e-2-''(--*)dx. (1) 
-OO 

This is simply the regular STFT with a modulation factor e^'^"'*, which will be convenient for our purposes. 
We will also occasionally write /i ^ /2 if the inequality Ci/i < /2 < (^2/1 holds, where Ci and C2 are 
constants independent of fi and /2. Finally, we will let sinc(x) :~ ^'"^^^-^ . 

We now consider a function / having the AM-FM form f{t) ~ A{t) cos{2TT(f>{t)). We want to deter- 
mine (/)'(t), which intuitively describes the local rate of oscillation of / at t. This leads to the following 
general definition. 

Definition 2.1. Suppose / is a superposition of K AM-FM components, having the form 

K 

/(t) = ^ylfe(t)cos(2^0fe(i)) (2) 
fe=i 

with (t>'i;{t) > for all k. Then the ideal instantaneous frequencies IIF(/, Ak, 4>k) are defined to be the set 
of functions {0^(i)}i<fc<K- 

Note that this concept only makes sense if all the component functions Ak and (jjk are known, and it is 
impossible to define the IF of an arbitrary function / in this way. In fact, an arbitrary / will generally not 
have a unique representation of the form ([2]), with many different choices of Ak, 4)k and K possible. Even 
if we restrict K = 1, there is in general no way to separate the amplitude factor Ai from the frequency 
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factor cos(27r(/)i) without having some additional information on /. 

In data analysis applications, we typically only know the full signal / and we want to obtain an approxi- 
mation to the IIF. One of the most commonly used approaches for doing this can be described as follows 
For an appropriate /, we define the operator P'^f = (/ + iHf) /2, where Hf ~ (— isign(?7)/(?y))'' is 
the Hilbert transform of /. is known as the Riesz projection in mathematics and as the single- sideband 
modulation or the analytic signal in the engineering literature. We then consider the Hilbert transform IF 
given by 

The motivation for this concept is that the function / is assumed to have the form ([2)) with a sin- 
gle AM-FM component, f{t) ~ cos(27r0(t)), and under some conditions, lEnfit) is a good ap- 
proximation to (jj'it)- The Bedrosian theorem states that if supp(A) and supp(cos(27r0)) are disjoint, 
then P^f{t) = A{t) ■ P+ cos(27r(/i(t)). If we additionally assume that supp(exp(27ri(/))) C [0,oo), then 
P+ cos{2tt (t){t)) = e2'^*'^(*), and we have Wnfit) = ^Im (2TTi(j)'{t) + ^) = (j)'{t). 

These conditions on A and </> are fairly restrictive and can be hard to verify for real-world signals, partic- 
ularly the requirement that supp(exp(27ri(/))) C [0,oo). Even when they hold, the computation of IFh is 
often sensitive to noise and numerical roundoff errors due to the Hilbert transform computation and the 
possible cancellation of zeros in the numerator and denominator of In spite of this, IFh lias proven 
to give meaningful results for signals arising from a variety of applications, and is in widespread use in 
data analysis. We refer to the papers jlT] and |2] for more details on the physical interpretation of IFh- 
In Section [3l we will propose an alternative approach based on different ideas to approximate the IIF of 
a signal. We return to IFh in Section |4] and show how it can be computed from nonuniform samples of a 
bandlimited signal. 

3 Synchrosqueezing with the Short-Time Fourier Transform 

Synchrosqueezing is an approach originally introduced in the context of audio signal analysis in [H] and 
was recently developed further in [S]. Synchrosqueezing belongs to the family of time- frequency reassign- 
ment methods and is a nonlinear operator that "sharpens" the time-frequency plot of a signal's continuous 
wavelet transform so that it provides more useful information about the IF. In contrast to the reassignment 
methods discussed in [2], Synchrosqueezing is highly adaptive to the given signal and largely independent 
of the particular wavelet used, and also allows the signal to be reconstructed from the reassigned wavelet 
coefficients. We refer to the paper [5] for more details. 

In this paper, we take a slightly different approach to Synchrosqueezing than in [S], based on the modified 
short-time Fourier transform. We will develop our theory independently of the results in [5] and show that 
this new Synchrosqueezing transform provides a way to estimate the IIF of a given function from discrete, 
nonuniform samples of the function. We will make a series of definitions leading up to our main theorem. 
We first define the following class of functions. 

Definition 3.1. (Intrinsic Mode Functions) 

The space of Intrinsic Mode Functions (IMFs) of type B consists of functions / having the form 

f{t) = A(0e2^*'^(*) (4) 
such that for some fixed e <^1, A and ((> satisfy the following conditions: 
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We then consider the function class B^.d, defined as follows. 
Definition 3.2. (Superpositions of IMFs) 

The space B^^d of superpositions of IMFs consists of functions / liaving the form 

K 



for some K > and fk = Ake/^^''"' € Be such that the (j)k satisfy 

inf0'fe(t) -sup0;,_i(i) > d. 



t 



t 



Our main results in this paper will be stated for B^^d functions. Intuitively, functions in B^^d are composed 
of several oscillatory components with slowly time- varying amplitudes and IIFs (we call these components 
IMFs, following jSj), and the IIFs of any two consecutive components are separated by at least d. This 
function class is fairly restrictive, but it turns out to be a reasonable model of signals that arise in many 
applications, and the conditions on the IMFs will allow us to obtain accurate estimates of the IIF of /. 
Note, however, that Be,d is not a vector space. We next define a closely related notion. 

Definition 3.3. (Impulse Trains) 

Let r > and {a„} G l°° with ||{an}||;°° < T^- The impulse train class is defined by 

f2 iT + an+i-an)S{t-Tn-an)f{t): f e B.A ■ 

We can treat elements of the impulse train class 'J^ as tempered distributions. As T — > 0, / G ^ " 

converges weakly to / G Be,d with respect to Schwartz test functions, so D^'j'^"^ can be thought of as a 
sampled version of the continuous-domain B^^d class. In applications, we are given a collection of nonuni- 
formly spaced samples of the form {f{tn)}, t,i = Tn + Un, where {a„} represents a small perturbation from 
uniform samples {Tn} with sampling interval T. This can be converted to an element in the class V'^^'j"'"^ 
by multiplying f{tn) by a factor tn+i — tn = T + a„+i — a„. 

For a given window function g G S, we can apply the modified short-time Fourier transform ([IJ to any 
/ G with the tempered distribution / acting on g(- - t)e"27ri^( -t) ^ 5 When \Vgfit,r])\ > 0, we 

define the instantaneous frequency information uif{t,r]) by 

-/(t,^) = ^^M^. (5) 
27rzV,fit,rj) 

The idea with ^ is that dtVgf is a first approximation to the IF of /, and dividing by Vgf "removes the 
influence" of the window g and sharpens the resulting time-frequency plot. We can use this to consider the 
following operator, which will be our main computational tool in this paper. 



4 



Definition 3.4. (STFT Synchrosqueezing) 

For f eV 
defined by 



For / G the STFT Synchrosqueezing transform with resolution a > and threshold 7 > is 



wliere (t,^) G M x aN and |-| denotes tlie Lebesgue measure on 



r,: \^~u;f{t,rj)\ < |, \Vjit,r,)\ > 0<V<^ 



(6) 



8°"'^ f is a kind of highly concentrated version of w/ that "squeezes" the content of w/ closer to the IF curves 
in the time-frequency plane. We are finally in a position to define an alternative notion of instantaneous 
frequency, based on the concepts discussed above. 

Definition 3.5. (Synchrosqueezing-based Instantaneous Frequency) 

The Synchrosqueezing-based IF with resolution a and threshold 7 of / G B^^^d^ estimated from a correspond- 
ing / e V'^j''"\ is a set-valued function lFsf{t) : R ^ 2"^ given by 



IF. 



= {an : n e N, an) > o} 



As a motivating example of this concept, let f{t) = e^'^** G B^^d- The IIF is clearly the singleton {1}. For 
the window g{t) = e"'^*^, we can compute Vgf{t, rf) = g^'^**"'^'''"^'^, so ?y) = 1 for all {t, rf). This means 
that for 7 = 0, S°''''f{t,^) will be supported on {[ija, [^Iq;} for all t, which is close to {1} for small a. 

We will show in Theorem 13.61 that if / G -p^^^""^ is a sampled approximation of / and a and 7 are chosen 
appropriately, then it has the same IF5 set. 

We will use several auxiliary notations in the statement and proof of our main theorem. Let 

K 00 

/W=E E - Tn - aJ(T + a„+i - a„)Afc(i)e2-^'=(*) G2?Jf"^ (7) 

k—1 n— — 00 

We then define the following expressions: 



J^) := / \u\"\g^^^^Hu)\du 

I := [J In, Tn — [Tn, Tn -\- a„] if a„ > or \Tn + a„, Tn] if a„ < 



riGZ 
K 



E, := J2'\\^'k\\L^ih+7r\\Ak\\L^l2) 

k=l 

^ ( 1 



k=\ ^ 

Ez := sup 110', 11^^ 

l<k<K 

This leads to the following results. 

Theorem 3.6. Let < T < 1. Suppose we have / G 'D^d'^"^ as in (O with ||{a,i}||;oe < T^, and a window 
function g G 5 with supp(g) C [— and Jj-\g{x)\ + \g'{x)\dx < y ||{an}||;oo for a constant k. Then 
there exist numbers E2 = E2{Ak, 4>k, g) and E'2 = E2{Ak, (f>k, 9) such that if we have a given resolution a 
satisfying 

2(E'. + E'n) 

- E1+E2 

then the following statements hold. 
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1. Let < r/ < ^ and fix fc, 1 < < K. For each pair {t,'q) e Zk := {{t,r]) : \r] - < |} witlr 
\Vgfit,rj)\ > E1+E2, we have |w/(t, ry)-0;(t)| < f . If (t,??) for any fc, then \Vgf{t,Tj)\ < E1+E2. 

2. Suppose we have a threshold 7 such that Ei + E2 < 7 < \Vgf{t, r/)| for all {t, if) G 2'^. Then for all i, 
S"''^f{t,^) is supported in the 2A'-point set Ui<fc<x{L^^Ja. T^^l"}- 

Theorem 13.61 savs that the IF5 defined by STFT Synchrosqueezing approximates the IIF set accurately up 
to the preassigned resolution a, without knowing anything about the symbolic form of /. The result also 
does not depend on the precise shape of the window g that we use, so the IF5 is in a sense adaptive to 
the structure of the sampled function /. The procedure suggested by Theorem 13.61 can be implemented as 
follows. We first discretize / on a uniform grid (finer than T) by zero-padding in between the impulses. 
For each t, we can then compute Vgf and dtVgf = —Vg'f + 2TTiriVgf on using FFTs and use the results 
to approximate 3°''"' f(t,an) for n G N. In practice, the upper bound on 77 in ([H]) is determined by how 
finely / is discretized, and as long as the threshold 7 is not too small, we ignore the locations where the 
denominator in ([5]) is close to zero and avoid numerical stability issues. We finally find the (numerical) 
support of S"''^ f{t, •) to determine the component(s) of lFgf{t). 

There is a tradeoff between a and 7 and the fluctuation of the IIF components, and this is a kind of uncer- 
tainty principle that is inherent to the IF5 concept. The IF5 is only meaningful up to the resolution a, and 
beyond that, we cannot approximate the IIF to any further level of accuracy. In the simplest case when the 
amplitudes {Ak} and IIF components {4>'k} ^-U constant, the lower bound on a will be small and will 
allow us to specify a fine resolution a, resulting in a very accurate estimate of the IIF. On the other hand, if 
the {A'f^} or {(t)'k} are large, then the lower bound on 7 will be significant, making S'^''^ f{t,(,) = for most 
t. In physical terms, the existence of a 7 in Theorem 13.61 ensures that the magnitude of each component is 
not too small and its IIF is not too rapidly oscillating, or else the component would be indistinguishable 
from noise. 

The proof of Theorem 13.61 involves a series of estimates. Let / and / be given as in Theorem 13.61 In what 
follows, we let Qk{t,r]) = Ak{t)e'^^^'^>'^*^ g{r] — 4''k{'t)) to simplify some notation. 

Lemma 3.7. If {t, rj) G Zk for some fixed fc, 1 < fc < A', then 

\Vgf{t,ij)-Qk{t,ri)\<E^ (8) 

and 

< E[. (9) 



If (i, 77) ^ Zk for any fc, then 

\Vgf{t,7^)\<E^ and 

Proof. Note that for any the bandwidth condition on g shows that Qz(t, rj) = whenever {t, if) ^ Zi. We 
assume (i, 77) G Zk for some fc, as the other case can be done in exactly the same way. For the first bound 
(O, Taylor expansions show that 

It follows that 

\Vg.f{t,ii)-Qk{t,rj)\ 



dtVgf{t,rj)-ct>'k{t)Qk{t,v) 



dtVgf{t,r,) 



<E[. 
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K 



k=l 



< 



K 

E 



k=l 



dx 



+\Mt)\ 



g{x - i)e-2"')(^-*) 



dx 



< 



e(/ k 

fc=l 



'i,\\^^\x~t\\g{x-t)\dx+\Au{t)\ / ^\\(l)l\\^^\x~t\^\g{x--t)\dx 



K 



< ^6||<^L.|Lo.(/l+7r||Afc||^^/2). 



k=l 

We follow similar arguments for the second bound We first have the estimate 

< \Ak{x)<p'^{x) - Au{t)<P't,{t)\ + |1 -e2-(-^^(-)+*''(*)+*^(*)(--*))||Afe(O0'feWI 

< (II^^IIlo. UIWl^ + PfclLoe ll^'.ll^^ \x-t{^ 

Thus we obtain 

1 



27ri 

r 

E 



dtVg.f{t,T^)-<P'k{t)Qk[t,r^) 



k=l 



K 

^ E 

k=l 



1 



) 

OO 
^ — OO 



1 
2^ 



K 



k=l 



Lemma 3.8. Let < ?/ < y. Then 



Vgf{t,,j)-Vgf{t,v) 



and 



1 

2^ 



< E'. 



□ 

(10) 
(11) 



for some numbers E2 = E2{Ak,4'k,g) and £"2 = E2{Ak,4>k, g) independent of rj or T. 
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Proof. Suppose h ^ C^. Then denoting ~ Tn + a^, we have 

J2 - tn)h{t„) - / h{t)dt 

J —OO 



E 



n— — oo 



h' {u){u — tn+i)du 



(12) 



/OO QO 
-OO „_ ^ 



n— — OO 



Shice Vgf{t,ff) — X!^-oo(^^^+l ~ tn)f{tn)9{tn — , we use the above calculation to find that 

\Vgf{t,rj)~Vgf{t,lj) 

OO 

n— — OO 

2er||0;,ll^^/o + 2rpfc||^^^V 



K 

^ E 

/c=l 

a: 

^ E 



G(u)C(M)du 



(13) 



where we let G{u) 2^i^fe(M)(0',.(«) -77)5(M-i)e'"(*''(")-"'' and C(u) Er=-oo("- Wi)X[t„,t„+i](«)- 
For brevity, we will fix k and omit the subscripts on and (j)k in what follows. The function C{u) is well 
approximated by the uniform sawtooth function 

TTT/ N / rr,/ / N >^ slu ( 27rnu/r) T 

Wiu) = ^ {u~Tin+l))x[Tn,T(n+i)]iu)^-22 Vn/T 2' 

n— — OO n— 1 



where the last equality holds for almost all u. We define the difference D{u) = W{u) — C{u). Since G ^ S, 
using the Parseval theorem gives 



G{u)C{u)du 



G{u) '^^^^^""^^^ dit - G(0)^ + / G(u)i?(7.)du 



< 



^ nn/T 



^|G(0)| 



G{u)D{u)du 



(14) 



We now consider each term in ([H)) separately. For the first term, we will need an estimate of G{C)- We 
claim that G will be insignificant outside intervals centered at (^'{u) — rj. For any J > 0, let ^ be such that 
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\4>'{u) — ?y — ^1 > J. Then we have 



\G{0\ = 2n 
= 2n 
< 27r 



A(u)((l)'{u) - 77)e2"('^(")-(''+«)")g(u - t)du 



2M'l>H-{n+i)u)Q [ Mu){(l}'iu) -r])g{u-t) 



J(e Pll^^ /o + e ll^'ll^^ /o + r,) + e ||0'||^^ \\A\\^^ h 



J2 



+?7 



,/2 



(15) 



where the constant Ci = Ci{Ak^ (pk^g) is mdependent of r/ or J. We take J ~ niax(^^ — i+, L — y) 
with the numbers L+ and chosen such that < \4>'{u) — 7?| < for all u. This gives 



E 



G 



2T 



47rnniax(^ - L+, L- - f ) 



< Ci(77 + l)f;2max(|L+-L-|,l)(^^^ 



Cimax(|L+ 



12 



(16) 



For the second term in (|14p . we consider two cases. If |0'(u) — 77] > J, then calculations similar to those in 
(fl31) show that 



|G(0)| 



27r 



(Mu)-nu)Q ( A'{^)9{u - + A{u)g'{u - t) 



J2 



1 1 



for some constant C2. On the other hand, if \4>'{u) ~ri\ < J, then we can easily estimate |G(0)| < ||G|[^i < 
27r JIq- We now choose a J that minimizes max(C2 (7 + jr) , 2tt \\A\\j^ac JIo) over all < 77 < y. 

This results in a bound of the form 



j\Gm<C3T 



(17) 



for a constant C3. The third term in p4p is controlled by the nonuniform perturbation {«„}. Note that if 
u^l, then |-D(m)| < 2 ||{an}||;oo, and for u G X, we have |-D(u)| < T + 2 ||{an}||;=o. This gives 



9 



G{u)D{u)du 

OO 

G{u)D{u)du+ I G{u)D{u)du 



< ll'^'-^llL=e||A||^^ (r+||{a„}||,^) / \giu)\du + Io\\{a,,}l 



< 



+ v) Uh^ ((T + T')kT + /oT^) . (18) 
Combining ([13]), (fT6|) . (|17|1 and p8)l . we end up with an estimate of the form 

Vj{t, v) - Vgfit, r;)| < C,{T + T\tj + 1)), 

for some C4 = C4{Ak, (l>k, g) independent of T, {a„} or 77. This finally implies the result ([T0|) . For the 
derivative inequality (|11|1 . we can simply use the result of ()10[) with g replaced by g' , which shows that for 
some C5, 



1 



2tt ' 

K 



dtVjit,f^)-dtVgfit,r^) 

duiAk{u)g{u - t))e2"('^^(")-''("-*)))C(M)du 



1 

2^ 



□ 



Proof of Theorem VJ.bX For the first part of Theorem l3.61 we suppose that (t, ?y) G for some fc, \ Vgf{t, if)\ > 



El + E2 and < < A. Putting Lemmas 13.71 and 13.81 together, we get Vgf{t,r]) — Qk{t,ri) 



and 



^A^gfit, v) - (t>'k{t)Qkit, v) <E[+ E!,. We conclude that 

^^dtVgf{t,rj) - , Mt)(Qk{t,v) - Vgf{t,Tj)) 



< 



< 



E[ + E'^ + EajEi + E2) 
El + E2 



Vgfit,v) 



< El + E2 



(19) 



For the second part of the theorem, let t be fixed and suppose ^ e aN, ^ ^ Ui<fc<if{ L"^ir^J"^7 T^TT^IQ!}- 
If {t,ri) ^ Zfc for any k, then by Lemmas 13.71 and 13.81 we would have \ Vgf{t,r])\ < Ei + E2 < j- But since 



-0'fc(t)- 



l^/(^j'7)l ^ 7 by assumption, (t, r/) G Zk for some /c. If ^ > 



a, then the estimate (fT9|) gives 
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and similarly, if ^ < [^^^Ja, we get a;/(t, 77)! > \iof{t,T])\—^> ^. This means that {77 : — a;/(t, 77)! < 
f , \Vgf{t,Tj)\ > 7, < 7? < i} = and so S"-y{t,^) =0. □ 

Remark. We have studied functions / G on the entire real line here, but we can instead consider such 
functions supported only on a finite interval, and the same results will hold with largely only notational 
changes. We can also consider real- valued functions / that have the exponentials e'^'^^'P^it) j-gplaced 
by cos(27r0fc(t)). and similar results will hold. The assumption that T < 1 in Theorem 13.61 was made for 
convenience in the proof, but it poses no loss of generality since lower sampling rates are equivalent to 
simply rescaling / (and thus also its IIF). 



Remark. In Theorem 13.61 '^e required the window g to be bandlimited to [— f , |]. This assumption is not 
strictly necessary and was made here to simplify the presentation. In general, it is enough for g € S to 
be such that \g\ is small outside [— |, |]. The Gaussian window g{t) ~ e~^^ 1^ works well in practice. 
For such a window, there will be an extra term in the error estimates in Lemma |3.7[ but the results still 
remain essentially the same. (e.g. if we assume ||g||^oo(jj^[_^ < e, then instead of having Qi{t,ri) = 
for {t,7]) ^ Zi, we would get \Qi{t,7])\ < Ce for a constant C.) 

We finally show that if a and 7 are chosen properly, the calculation of IF5 is robust to sample noise. The 
argument is almost the same as in the proof of Theorem 13.61 



Theorem 3.9. Let < T < 1. Suppose / € T^f'j'^"^ is as in ^ and let g be a window with the same 
conditions as in Theorem l3.6l Assume that the samples {f{tn)} = {f{Tn + an)} are contaminated by noise 
{Nn} with ||{iV„}||,oo < T, i.e. we are given / = / + I]^^_oo(Wi - tn)d{- - tn)Nn. Suppose we have a 
resolution a satisfying 

^ 2{E[ +E',+Io + 21', + ^(/^ + 2/^)) ^ ^ ^ 
E,+E,+I, + 2I', +2^^' 

where E2 and E'2 are as in Lemma |3^ Then the following statements hold. 

1. Let < 77 < ^ and fix k, l<k< K. For each pair (t, 77) e Zu with \Vgf{t, ri)\ > Ei + E2 + Iq + 21,, 
we have |w/(i, 77) - ^'^(t)! < f . If (t, 7/) ^ Zu for any fc, then \Vj{t, 77)] <E^+E2 + h + 21',. 

2. Suppose we have a threshold 7 such that Ei + E2 + lo + 21, < 7 < \Vgf{t,ri)\ for all (t, 7/) 6 Zk- 
Then for all t, S"''^f{t,^) is supported in the 2A'-point set Ui<fc<x{L^^Ja: [^^1"}- 

The only part of the proof different from Theorem 13.61 is the estimate in Lemma 13. 8| which we replace by 
the following inequalities. 

Lemma 3.10. Let < 77 < i. Then 

\Vgf{t, v) - Vgfit, v)\<E2 + Io + 21', 

and 

^ \dtVgJ{t,r^) - dtVgf{t,7^)\ <E'2+Io + 21', + ^{I', + 21'^). 
Proof. Lemma 13.81 and the inequality (|12p imply that 

\VgJ{t,ri)~Vgf{t,l^)\ 

< Vg~f{t,r^)~Vgf{t,r^)\ + \\N^l^ ^ |(t„+i - i„)ff(i„ - t)| 

< E2+T{h + 2TI',), 



n— — QO 
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and similarly, 



— \dtVj{t,ij)-dtVgfit,Tj)\ < i?^+7?T(/o + 2r/^) + — T(/^ + 2T/^')- 

□ 

Combining these bounds with the other estimates in the proof of Theorem 13.61 gives Theorem 13.91 

4 A Bandlimited Reconstruction Approach 

We now consider a somewhat different formulation and approach towards our problem, based on the 
traditional IFh concept. Suppose we have the samples {f{tk)} of a function / G with supp(/) C [—b, b]. 
We want to determine / and use it to find IFnf- In the case of uniformly spaced samples, i.e. tk — Tk for 
a constant T > 0, it is well known that / can be recovered if the sampling rate i exceeds 26, the Nyquist 
frequency of /. There are analogous results for very general classes of nonuniform sampling points {tk\ 
[T1I12]. As in Section [HI we will consider sampling points {tk] ~ {Tk + ak] that are small perturbations of 
uniformly spaced points. Suppose the sequence {tk\ is indexed so that < tfe+i and that for some T > 0, 

sup|ifc-Tfc| < ^. (20) 
fc 2 

If ^ > 26, it can be shown that Ci ||/||^2 < ||/(ifc)||/2 < (^2 ||/||2,2 , where the constants Ci and C2 depend 
only on {t^} and 6 [8l[T0]. When the tk are uniformly spaced, i.e. tk = Tk, the Plancherel formula shows 
that Ci — C2 = r~^/^. For the rest of this section, we assume that {tk} satisfies ((20|) with ^ > 26, 
which intuitively means that the sampling points {tk} are evenly spread out over R and have an "average" 
sampling rate above the Nyquist rate. 

One of the standard approaches for recovering / from the samples f{tk) involves solving the linear least- 
squares problem 

2 

(21) 

P 

where {h(n,t)}n^[-N,N] is a given set of basis functions and the Wk are some weights such that [wkl ~ 1 
for all k \10 \ \T2 \ \9\. If the basis and weights are chosen appropriately, the function 

N 

fN{t) = CnMh{n,t) (22) 

is a good approximation to frit) f{Tt). There have been a couple of efficient algorithms developed 
around this idea, usually using the DFT basis h{n,t) = {2N + i)-i/2g-27r«it/(2JV+i) ^^le weights 
Wk ~ tk+i-tk-i ^ -pQY these choices, as N ^ cxd, /n converges uniformly to /r on compact subsets 
of M, and the matrix computations used in solving the problem ([?!]) are well-conditioned and can be 
accelerated using FFTs [10]. Once we have a expansion of the form ([22| . it follows by linearity that 
P+/7v(i) = J2n=-N ^n,NP^h{n,t), and we can use this to find IFnf- 

We will prefer to use the basis h{n,t) = smc{t — n — M) in this paper, where M is an integer constant, 
along with the same weights Wk as above. In this case, /n converges to fx uniformly on the entire 
real line, and we have found that in practice, this results in a better accuracy with the computation of 



N 



{ck,N} = argmin 

{cfc}eC2«+i 



Wk 



.f{tk)- ^ Cnh{n,tk/T) 



n=-N 
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P^fpf. This choice of h{n,t) can be justified by the fohowing argument. For any vector {c^-} £ with 
supp(cfc) C [M - N,M + N], we have 

TV 

hit) - ^ c„+Msinc(< -n- M) 



ll/T(fc)-Cfe||p 



Ti=-TV 
TV 



71= -N 



TV 



f{tk)- J2 Cn+MSmc{tk/T-n-M) 



n=-N 



By taking the minimum over {cfc}, it follows that for {ck.N} given by (|2ip . 

Wfrik) -Cfc,Ar||;2 ~ ||/T(fc) - /T(fc)X[T\/-TV,A/+TV](fc)| 

as oo. We can then conclude that 



-> 



II/t-/ 



TV 



< 



fT~L 



TV 



,^ < ||/T-/Tv||i. = WMk) - dk.NWp 



So we essentially determine uniform samples of / from the nonuniform samples {/(tfc)} by solving the 
problem pip , and use those in a classical sampling series to compute /. In practice, we consider a finite 
number of samples spread over some interval [Ji, J2], and we set AI = [-^{Ji + J2)J to center the basis 
h{n,t) appropriately. 



Once we have recovered the function /, we perform some elementary calculations with Fourier transforms 
to find that 

TV 



-P^/Tv(i) = ^ c„,Arisinc 



t-n- M 
2 



71= -N 

and 

(j + ^(i_„_M))e-(*-"+*^) 
-P hit) ™ 2. ^".A' 2n{t-n-Mr ' 

n=-N ^ ' 

from which we can approximate IF///. Note that the Synchrosqueezing-based method discussed in Section 
[3] determines the IF components (as defined by IF5) of the function / directly, while the approach consid- 
ered in this section treats the function / as a whole and determines its IF (given by IF/f ) after the signal 
itself has been recovered. 



We finally make a few comments on the applicability of the framework discussed in this section. The 
assumption that f € can be weakened to / G L°° in practice. Suppose _B is a function such that B is 
smooth, B{0) = 1 and supp(i?) C [—1, 1]. For y > 26, any / e L°° with supp(/) C [—6, b] can be expressed 
as the sampling series 

f{t)^ /(Tfc)sinc(^|-fc)i?(^(^l-2fe)(t-fc)), 

which converges uniformly on compact sets If we restrict our attention to a fixed finite interval /. 
then / can be well approximated on / by taking a finite part of this series fp S L^, which has supp(/p) C 
[— (^ — &), (^ — 6)]. As long as / is oversampled, our reconstruction method will approximately recover 
fp. In the next section, we will apply the least-squares method to AM-FM functions / of the type in 
Definition 13.21 which are generally not strictly bandlimited, but by considering Fourier series on a given 
finite interval /, such functions can be closely approximated on I by bandlimited L°° functions that have 
a sufficiently high bandwidth (see [14] for details). In practice, if our sampling interval T is greater than 
supfc IIF(/, Ak, (pk), we can recover / effectively. 
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5 Numerical Experiments and Applications 



We use the algorithms discussed in Sections [3] and 2] on several test cases. We consider an AM-FM signal, 
a chirp signal, a bandlimited Bessel function, a Fourier harmonic contaminated by noise, an undersampled 
harmonic, a multi-component signal and finally, a signal with interlacing IIF elements. These are respec- 
tively shown in Figures [TJI] below. We compute the IF5 of these signals using STFT Synchrosqueezing and 
the IFh using the bandlimited reconstruction method. For the IF5 computation, the Gaussian window 
function g{t) = q-^-^"^*- and the parameters a = 0.1 and 7 = 10~^ work well in practice, and we use these 
values unless stated otherwise. 



In all of these examples, we use sampling times that are random perturbations of uniformly spaced times 
and have the form i„ ~ Tn -\- r'a„, where T' < T and {a„} is a fixed realization of a white noise pro- 
cess uniformly distributed in [0, 1]. For purposes of comparison, we also test most of these examples with 
uniform samples as well, i.e. T' = 0. In the figures below, the first two images are respectively IF5/ and 
IF/// determined from uniform samples, while the two latter images are IF5/ and IF^f for the general 
case of nonuniform samples. 




Figure 1: The AM-FM signal f{t) = (2 + cost) cos(27r(3i + cost)), t e [0, 30]. The samples are taken at 
tn = O.ln + T'an, T' ^0 (left images) and T' = 0.08 (right images). The IIF is 3-sint, and both methods 
produce results that are reasonably close to this. 




Figure 2: The chirp signal f(t) = cos(27r(t + 0.05t^)), t € [0, 30], with samples taken at t„ = O.ln + T'a„, 
T' — (left images) and T' = 0.08 (right images). The IIF in this case is 1 -I- O.lt, and both methods 
produce results that are very close to this. 
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Figure 3: The bandlimited signal f{t) = Jo(67r(f; — 10)), t G [0,30], where Jq is the Bessel function of 
order (see [15]). The samples are taken at t„ = O.ln + r'a„, T' = (left images) and T' = 0.08 (right 
images). There is no way to determine the IIF in this case, but it can be shown that Jo(0) = 1 and 

Jo{t) ~ \l cos(|t| — ^) when \t\ is large, so we would expect the IF to be roughly the constant 3. The 

computed IF5/ agrees with our intuition here. On the other hand, IF/// is highly oscillating, particularly 
around t = 10, which indicates that it is unable to clearly distinguish between the amplitude and frequency 
factors of /. We also show the graph of / itself at the bottom for clarity. 




Figure 4: The signal f{t) = cos(87rt) + Nt, t £ [0,30], where Nt is a realization of Gaussian white noise 
with mean and variance = 0.4. The samples are taken at i„ = O.ln + r'a„, T' = (left images) and 
T' — 0.08 (right images). The computation of IF^/ is fairly robust to the noise, while that of IF/// gives 
very poor results. 




Figure 5: The signal f{t) — cos(107r<), t e [0, 30], with samples taken at t„ — 0.25n + 0.2a„. The IIF is the 
constant 5. Note that this signal is heavily undersampled, with a sampling rate of less than half its Nyquist 
rate. We take 7 = 6 for the STFT Synchrosqueezing computation, which produces a good result for IF5/ 
despite the low sampling rate, but the bandlimited reconstruction method cannot determine IF///. In 
this example, we also show the time-frequency plot of 3°'''^ (the first image) to illustrate how the IF curve 
appears in it, before we extract it out. 



15 




Figure 6: The two-component signal f{t) = cos(27r(2i + 0.2 cost)) + cos(27r(3t + OMt"^)), t £ [0, 30], with 
samples taken at t„ = O.ln + T'an, T' — (top images) and T' = 0.08 (bottom images). We would like 
to recover both elements of the IIF set {2 — 0.2 sint, 3 + 0.04t}, and the IPs calculation succeeds in doing 
this. In contrast, the IFh concept cannot separate the components, instead roughly giving their average 
2.5 + 0.02t — 0.1 sint, and its computation also exhibits spurious singularities. 




Figure 7: The signal f{t) ^ cos(57rt)+cos(27r(i-|-0.05t^)), t G [0, 30], with samples taken at t„ = 0.l77,+T'a„, 
T' = (top images) and T' = 0.08 (bottom images). This example is similar to the previous one but the 
IIF curves cross each other. However, the plot of S"''^ shows a sharp separation between the two curves 
around the crossover point and can distinguish between them clearly. 

The examples in Figures [3]l7] show the advantages of the IF5 approach over the traditional Wh concept. In 
Figure[5l the bandlimited reconstruction method fails to recover / accurately due to the low sampling rate, 
and this is reflected in the computation of IF///, but STFT Synchrosqueezing still manages a good result. 
This can be possibly explained by the fact that the Nyquist rate is essentially a concept for bandlimited 
signals and is only relevant in that context, whereas our STFT Synchrosqueezing theory is built around 
Be^d signals and only estimates their IIF, not the signals themselves. In Figures IH [6] and [71 our sampling 
rate is high enough and the bandlimited reconstruction method can accurately determine / itself, but the 
subsequent calculation of IF/// amplifies the effects of any noise or numerical roundoff errors. In contrast, 
we find that IF5/ is robust to such disturbances. We also note that determining a meaningful IF for the 
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type of signal in Figure [7] is often difficult [7] and such signals are certainly not in the class Be^d, but the 
result is nevertheless very good. 



We now discuss a real-world problem in electrocardiography (ECG) to which our methods are applicable. 
In addition to describing the heart's electrical activity, the ECG signal contains information about a signal 
describing respiration. It is important in many clinical situations to be able to determine properties of this 
respiration signal from the ECG signal. For example, in an examination for tachycardia during sleeping, 
where only the ECG signal and no respiration signal is recorded, it allows for the detection and classifica- 
tion of sleep apnea. It is well known in the ECG field that the customary surface ECG signal is influenced 
by respiration, since inhalation and exhalation change the thoracic electrical impedance, which suggests 
that the respiration signal can be estimated from the ECG signal. The ECG-Derived Respiration (EDR) 
class of techniques [13], in development since the late 1980s, accomplish this and have proved to be a useful 
clinical tool. We now show that STFT Synchrosqueezing provides an alternative method to extract key 
information about the respiration signal from an ECG signal. We can find the instantaneous frequency 
profile of the respiration signal, which gives a more precise and adaptive description of respiration than 
many of the existing techniques. 

In Figure [8l we are given the lead II ECG signal and the true respiration signal of a healthy 30 year 
old male, recorded over an 8 minute interval. The sampling rates of the ECG and respiration signals are 
respectively 512Hz and 64IIz. We take the R peaks (the sharp, tall spikes in the first image in Figure [5]) 
of the ECG signal and use these samples to approximate the IF of the respiration signal, without using 
any knowledge of the actual respiration signal. We do not have samples of the respiration signal itself, 
but we can view the R peaks as samples of an envelope of the ECG signal, and based on the physiological 
facts discussed above, this envelope would be expected to have the same IF profile as the actual respiration 
signal. We apply the methods from Sections |3] and U] to compute IFs and IFh from the R peaks, and 
use the actual, recorded respiration signal to compare the validity of our results. Let the ECG and true 
respiration signals be denoted by E{t) and R{t) respectively, where t G [0,480] seconds. There are 589 R 
peaks appearing at times tk G [0,480], tk < tk+i, 1 < fc < 589. For the calculation of IF5, we use the 
impulse train-like function 



and for IFh, we simply use the samples {E{tk)}. The results are shown in Figure |8] below. The third image 
in Figure [8] shows that the IF5 computed from fRpeaks is a good approximation to the IF5 of the true 
respiration signal R{t). On the other hand, the IFh determined from the R peaks in the fourth image has 
little in common with IFnRit). In fact, IF//i?(<) is often negative and admits no obvious interpretation, 
unlike the IFs- This is likely a result of the fact that ECG measurements usually contain large amounts of 
noise, which IFh does not handle well, and many standard noise reduction techniques are not applicable 
here since they would smooth out the R peaks. In Figure [^1 it can be seen that the spacing of respiration 
cycles in R{t) is refiected by IF5 of the R peaks; closer spacing corresponds to higher IF5 values, and wider 
spacing to lower IF5 values. 





otherwise 
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Figure 8: First Image: A 10-second part of the ECG signal. The R peaks are highlighted in red. Second 
Image: A 10-second part of the respiration signal. Third Image: The IF5 computed from fnpeaks (blue) 
and the IF5 of the actual respiration signal R{t) (red). Fourth Image: The IFh computed from E{tk) 
using bandlimited reconstruction (blue) and the IF/f of R{t) (red). 
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Figure 9: The first 300 seconds of R{t) (blue) with the IF5 estimated from fupeaks (red) superimposed 
on top of the graph of R{t). 
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